#! /bin/sh

# define parent directory
MYHOME=.

# output files
wfp=./pressure
wfr=./radial
wfz=./vertical
outf=./info

# input parameters
tsec=2.048		# length of computed trace (seconds)
fs=0.07			# sampling paramter 0.07<fs<0.12
decay=50.0		# decay factor to avoid time series wraparound
flt=0			# earth flattening flag
lsource=3		# layer on top of which the source is located
nlayers=6		# number of layers to process
nw=60			# number of frequencies to process for ref. modeling
nor=1			# number of receivers
dt=0.004		# time sampling interval (secs)
p2w=3.8			# maximum ray parameter (s/km)
lobs=3			# layers on top of which the receivers are located
bx=0.025		# beginning range (in km)
fx=2.975		# final range (offset) (in km)
dx=0.05			# range increment (in km)
pw1=0.0			# begin of low end Hanning window for ray parameter s/km
pw2=0.0			# end of low end Hanning window for ray parameter s/km
pw3=0.0			# begin of high end Hanning window for ray parameter 
pw4=0.0			# end of high end Hanning window for ray parameter s/km

#source parameters
stype=2			# source type.=1 for explosive, 2 for fault plane
h1=1.0			# vertical linear part of the source
h2=0.0			# horizontal linear part of the source
m0=0.0			# seismic moment
m1=0.0			# [1][1] component of the moment tensor (stpe=1)
m2=0.0			# [1][2]=[2][1] component of the moment tensor (stype=1)
m3=0.0			# [2][2] component of the moment tensor (stype=1)
delta=0.0		# dip in degrees (stype=2)
lambda=0.0		# rake in degrees (stype=2)
phis=0.0		# fault plane azimuth in degrees (stype=2)
phi=0.0			# receiver location azimuth in degrees (stype=2)

#parameters needed only in layer interpolation is required
nlint=0			# number of times layer interpolation is required
nintlayers=		# array[nlint] of number of layers to interpolate
intlayers=		# array[nlint] of layers on top of which to interpolate
intlayth=		# array[nlint] of layer thicknesses to interpolate

# parameters required only if random velocity layers are requested
rand=0			# =1 to request random velocity layers
nrand_layers=0		# maximum number of random layers allowed
layer=0			# layer on top of which the random velocity layers
			# are inserted
zlayer=0		# thickness of random layers
sdcl=0			# standard deviation for compressional velocities
sdct=0			# standard deviation for shear velocities

# parameters required only if qoption is requested
qopt=0			# =1 to request the q-option
layern=0		# layer on top of which the q-option is invoked	
wrefp=1.0		# reference frequency for compressional velocities
wrefs=1.0		# reference frequency for shear velocities
epsp=0.001		# reference amplitude for compressional velocities
epss=0.001		# reference amplitude for shear velocities
sigp=0.1		# xxxx for comporessional velocities
sigs=0.1		# xxxx for shear velocities

# key input parameters, (can be input from files, clfile, ctfile, etc) 
cl=0.400,0.600,0.800,1.200,1.500,3.200
ct=0.300,0.400,0.600,0.800,1.000,2.500
ql=15.0,25.0,30.0,50.0,100.0,200.0
qt=7.5,12.5,15.0,25.0,50.0,100.0
rho=1.800,1.800,1.900,2.100,2.300,2.500
t=0.003,0.005,0.010,0.020,0.050,1.00

# output parameters
wtype=1			# =1 for PSV, =2 for SH
wfield=2		# =1 for displacement, 2 for velocity, 3 for accelera
vsp=0			# =1 for vsp, 0 for surface seismogram
win=0			# =1 for windowed frequencies
wavelet_type=2		# =1 for spike, =2 for ricker1, =3 for ricker2,
			# =4 for an akb wavelet
verbose=3		# =0 for no processing information, =1 for processing
			# information in the screen, =2 for processing 
			# information to a file (outf) =3 for both
red_vel=0		# reducing velocity, equal to maxcimum cl if set to 0
fpeak=25		# peak frequency for ricker or akb wavelet
nf=256			# number of frequencies in the output data
nt=512			# number of time samples in output traces
tlag=0.0		# time lag in seconds
nx=60			# number of output traces


# output filter parameters (can be input from files, filtypefile,dbpofile, etc)
nfilters=2		# number of filters to apply to output traces
filters_type=2,1	# =1 for high cut filter, =2 for low-cut filter
			# =3 for notch filter
filters_phase=0,0	# =0 for zero phase filter, =1 for minimum phase filter
dbpo=24,24		# DB/octave for filter slopes
f1=3.0,30.0		# frequency to start filter action
f2=0.0,0.0		# high end frequency filter

# plotting parameters:
agc=0			# apply agc to the plot 
wagc=0.5		# agc window length (sec) 
pbal=1			# trace balancing by rms value
qbal=0			# trace balancing by clip value
qclip=1.0		# clip by quantile on absolute value of trace
perc=99			# perc to apply to plot


##############################################################################
sureflpsvsh fs=$fs decay=$decay flt=$flt lsource=$lsource nw=$nw nor=$nor \
	tsec=$tsec dt=$dt lobs=$lobs bx=$bx fx=$fx dx=$dx nx=$nx \
	pw1=$pw1 pw2=$pw2 pw3=$pw3 pw4=$pw4 p2w=$p2w \
	stype=$stype h1=$h1 h2=$h2 m0=$m0 m1=$m1 m2=$m2 m3=$m3 delta=$delta \
	lambda=$lambda phis=$phis phi=$phi \
	rand=$rand nrand_layers=$nrand_layers layer=$layer zlayer=$zlayer \
	sdcl=$sdcl sdct=$sdct \
	qopt=$qopt layern=$layern wrefp=$wrefp wrefs=$wrefs epsp=$epsp \
	epss=$epss sigp=$sigp sigs=$sigs \
	cl=$cl ct=$ct ql=$ql qt=$qt rho=$rho t=$t nlayers=$nlayers \
	wtype=$wtype wfield=$wfield red_vel=$red_vel nf=$nf win=$win vsp=$vsp \
	wavelet_type=$wavelet_type verbose=$verbose nfilters=$nfilters \
	filters_type=$filters_type dbpo=$dbpo f1=$f1 f2=$f2 nt=$nt tlag=$tlag \
	fpeak=$fpeak filters_phase=$filters_phase \
	wfp=$wfp wfr=$wfr wfz=$wfz outf=$outf

# plot the results
suxwigb < $wfp title="####PSV pressure field perc=$perc f=$fpeak" perc=$perc \
	label1="time (s)" label2="Trace Number" &

suxwigb  < $wfr title="####PSV radial field perc=$perc f=$fpeak" perc=$perc \
	label1="time (s)" xbox=600 label2="Trace Number" &

suxwigb < $wfz title="####PSV vertical field perc=$perc f=$fpeak" perc=$perc \
	label1="time (s)" xbox=300 ybox=300 label2="Trace Number" &

exit
